Scalar Field Dark Matter: non-spherical collapse and late time behavior 
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We show the evolution of non-spherically symmetric balls of a self-gravitating scalar field in the 
Newtonian regime or equivalently an ideal self-gravitating condensed Bose gas. In order to do so, 
we use a finite differencing approximation of the Shcrodinger-Poisson (SP) system of equations 
with axial symmetry in cylindrical coordinates. Our results indicate: 1) that spherically symmetric 
ground state equilibrium configurations are stable against non-spherical perturbations and 2) that 
such configurations of the SP system are late-time attractors for non-spherically symmetric initial 
profiles of the scalar field, which is a generalization of such behavior for spherically symmetric initial 
profiles. Our system and the boundary conditions used, work as a model of scalar field dark matter 
collapse after the turnaround point. In such case, we have found that the scalar field overdensities 
tolerate non-spherical contributions to the profile of the initial fluctuation. 

PACS numbers: 95.35.+d 98.62.Gq 05.30.Jp 04.40.-b 04.25.Dm, 



I. INTRODUCTION 

Recently, scalar fields have played different roles in sev- 
eral scenarios related to astrophysical phenomena. The 
reason is that such fields are quite common in theoretical 
physics, specially branches related to theories beyond the 
standard model of particles, high dimensional theories of 
gravity an tensor-scalar theories alternative to General 
Relativity. In the present research we deal with the scalar 
field dark matter model (SFDM), which assumes the dark 
matter to be a classical minimally coupled real scalar field 
determined by a cosh-like potential. Such potential be- 
haves exponentially at early stages of the universe and 
as a free field (quadratic potential) at late times, which 
provides the field with the necessary properties to mimic 
the behavior and successes of cold dark matter at cosmic 
scales. In fact in 0.013 it was shown that the mass pa- 
rameter of the scalar field gets fixed by a desired cut-off of 
the power spectrum, which has two effects: i) the theory 
gets fixed and ii) there is no overabundance of substruc- 
ture, which standard cold dark matter cannot achieve. 
One important consequence is that the boson has to be 
ultralight with masses around m ~ 10~ 21 '~ 23 eV. This 
is a substantially important bound, because in the stan- 
dard dark matter models there are no such ultralight dark 
matter candidates. The benefit obtained however, is two 
fold: the scalar field can represent a Bose Condensate 
of such ultralight particles and the Compton wavelength 
forbids the scalar field to form cuspy structures. In fact, 
in it was shown that in order for the scalar field in- 
teraction to become long range, the system only needs to 
be condensed. 

After the fluctuation analysis about this candidate 
and its corresponding concordance with observations - 
the model mimics the properties of the ACDM at cos- 
mic scales-, the next step has to be in the direction of 
the study of structure formation and the explanation of 
local phenomena, like rotation curves in galaxies. Fortu- 



nately there have been important advances in such direc- 
tion ja, @|- About the gravitational collapse, in 7] it was 
shown that relativistic self-gravitating scalar field config- 
urations can be formed when they have galactic masses 
provided the mass of the boson is ultralight. Neverthe- 
less, because the gravitational field in galaxies is weak, 
the race turned into the newtonian limit of the system 
of equations, which was developed in 8J. The price to 
be paid is that it is not possible to apply the approach 
at very early stages of the evolution of the universe, and 
the profit is that the scalar field in the non-relativistic 
regime provides a clear interpretation within the Bose 
Condensate formalism and classical Quantum Mechan- 
ics. In both cases, the strong gravity and the Newtonian 
regimes, a wide range of arbitrary spherically symmet- 
ric initial configurations collapse and form gravitationally 
bounded and virialized objects with a smooth density ev- 
erywhere called oscillatons (except those that are related 
to unstable initial configurations that collapse into black 
holes in the strong field regime) 0, ^| . This property 
seems to be fundamental in order to form galactic ha- 
los, because several high resolution observations are con- 
sistent with regular galactic dark matter profiles in the 
center of the galaxies ^lj> which implies that this type 
of condensate could be an alternative to solve the cuspy 
density profiles of dark halos. 

Two pieces of the model that are in progress are the 
high energy fully relativistic case, including exponential- 
like scalar field potentials and the free field weak energy 
newtonian, both are complementary and necessary to ex- 
plore the SFDM hypothesis. Inspired in a late-time as- 
trophysical scenario, at stages after the turnaround point 
where weak field applies and the field is free, the question 
is whether dark matter halos are gravitationally bounded 
objects of scalar field which have been formed through a 
gravitational collapse of initial scalar field overdensities. 
The newtonian version of the Einstein-Klein-Gordon sys- 
tem of equations is provided by the Schrodinger-Poisson 



2 



equations (SP), which are the ones that lead the gravita- 
tional collapse of the system. This approximation should 
work for the evolution of an initial density profile af- 
ter the epoch when the overdensity fluctuation starts to 
evolve independently of the cosmic expansion. 

In the recent past, it has been found that in spheri- 
cal symmetry the SP system has equilibrium solutions of 
two types: stable, for which the wave function is node- 
less (called sometimes ground state configurations) and 
others, for which the wave function has nodes (called 
sometimes excited configurations) that decay into ground 
state solutions. From these solutions only the ground 
state ones are stable 0, 0|; even further, it has been 
found that such configurations behave as late-time at- 
tractors for initially quite arbitrary spherically symmet- 
ric density profiles |sl.ll3|. 

In Q it was shown that free scalar field overdensities 
after the turnaround virialize and tend to form ground 
state configurations. Due to a very general scale in- 
variance of the Schrodinger-Poisson system of equations 
such results were also valid for a structure of arbitrary 
mass 0- In ^3 was shown that the spherical collapse 
of SFDM tolerates the introduction of a self-interaction 
term, which on the other hand is associated to the self- 
interaction term of a self-gravitating Bose-condensate 
[l4| . and helps at allowing diverse sizes and masses of 
the final configurations. 

In this paper we go a step forward and study the col- 
lapsing process of non-spherical initial configurations, in 
particular those involving at most quadrupolar terms. In 
order to achieve this goals an axisymmetric code that 
solves the SP system is needed. We choose to deal with 
cylindrical coordinates denoted by (x, z), where x is the 
radial coordinate and z the axial coordinate. The SP 
system for the free field case in these coordinates reads 
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where h = c = 1 and we are using the rescaled vari- 
ables x — > mi, z — > mz, t — > mt and the wave function 
ip — ► V 'AnGip. This set of coupled partial differential 
equations in two spatial dimensions plus time is the core 
of the present manuscript. As mentioned before, these 
equations have stable solutions in spherical symmetry 
(see [3 for solutions also including a non-linear term in 
the Schrodinger equation); such solutions have shown to 
be not only stable, but also late time attractors for quite 
arbitrary initial density profiles Il3j. Therefore, our 
main task in the present manuscript will be to show that 
such solutions are still attractors even for non spherically 
symmetric initial density profiles. 

In the next section we present the code we constructed 
for the present purpose. In section IIIII we show how 
a ground state configuration reacts under non-spherical 
perturbations. In section llVl we show the evolution and 



fate of non-spherical initial density profiles. Finally in 
section ]V\ we draw some conclusions. 



II. THE CODE 



A. Description 



The present code is built under the same numerical fi- 
nite differences method used in the spherically symmetric 
case in [9lll3|. We approximate the continuous equations 
(11121 using centered finite differencing for both coordi- 
nates x and zona uniform grid defined by x = pAx and 
z = qAz, p, q integers; we used the same resolution in 
both directions (Ax = Az). The spatial differential oper- 
ator is the same in both equations and we dealt with both 
in the same manner: aside of the usual finite differenc- 
ing expression for the space derivatives only two delicate 
items were included related to the first order derivative 
with respect to x in 11121) : i) we staggered the grid in the 
x-direction in order to avoid the divergence of such term 
and ii) we transformed such term into = 2-^%, with 

' x ox ax z 7 

the last expression a derivative with respect to x 2 . 

Schrodinger equation. In the present case this is the 
evolution equation of the system. We discretize time t = 
nAt, n an integer and At the resolution in time. We 
solve this equation using a second order accurate explicit 
time integrator, which is a modified version of the usual 
three steps iterative Crank-Nicholson method pj|. 

Instead of using a characteristic analysis of the propa- 
gating modes to set an open boundary at the edges of the 
domain, we decided to use a sponge in the outermost re- 
gion of the domain. The sponge is a concept used success- 
fully in the past when dealing with the Schrodinger equa- 
tion (for detailed analyses see 0, El). This technique 
consists in adding up to the potential in the Schrodinger 
equation an imaginary potential. The result is that in the 
region where this takes place there is a sink of particles, 
and therefore the density of probability approaching this 
region will be damped out, with which we get the effects 
of a physically open boundary. 

Poisson equation. Equation J2J) is an elliptic equa- 
tion for U which we solve using the 2D five-point sten- 
cil for the derivatives and a successive over-relaxation 
(SOR) iterative algorithm with optimal acceleration pa- 
rameter (see e.g. [13 for details about SOR). In order to 
impose boundary conditions we made sure the bound- 
aries were far enough for the mass M = J \ip\ 2 d 3 x to be 
the same along the three faces of the domain and used 
the monopolar term of the gravitational field; that is, we 
used the value U — —M/r along the boundaries with 
r = \/x 2 + z 2 for the gravitational potential. At the axis 
we demanded the gravitational potential to be symmetric 
with respect to the axis. 
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B. Construction of equilibrium initial data 

Initial data. The nature of the SP system allows one 
to have plenty of freedom about choosing the initial data, 
that is: once we choose an initial wave function ip we inte- 
grate Q at initial time; this means that the initial wave 
function is quite arbitrary. In fact we implement this type 
of initial data for the case of non-spherical collapse. Nev- 
ertheless, for the purpose of testing our numerical tech- 
niques, we decided to use initial data corresponding to 
spherically symmetric ground state configurations. Here 
we briefly describe how these data are obtained. 

In spherical symmetry equations read: 



A is a scaling parameter. Property I|7I8|) implies that if a 
solution is found for a given central field value 0(0) = 0o 
(e.g. (f)(0) = 1 p(0) = 1) it is possible to build the 
whole branch of ground state equilibrium configurations. 
For instance, if the plot M vs p(0) is to be constructed, 
we know from @ that for 0(0) = p(0) = 1 we have M = 
2.0622; using the relations I|7I8|I for and M we find 
A = (p/p) 1 ^ 4 — M/M, which implies that the desired 
plot is given by the function M = 2.0622 (p(0)) 1/4 for all 
central values of the scalar field density. This function 
is used later on when showing the attractor behavior of 
these configurations. 
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Uip 



dr 2 



(rU) = rifitp* 



(3) 
(4) 



where r = va 



is the spherical radial coordinate. It 



is assumed a time dependence of the type if) = (f>(r)e tu ' t , 
and demand the conditions of regularity at the origin 
(f)(0) = d x (j)(0) = and isolation (f>(x — » oo) = 0, the sys- 
tem becomes an eigenvalue problem where the frequency 
of the wave function is the eigenvalue. The system to be 
solved reads 



dr 2 



(rU) 



2r(U - w) 
r(f) 2 . 



(5) 
(6) 



We use a shooting method that bisects the value of u> for 
a given central value of (f> that satisfies the boundary con- 
ditions. That is, one constructs a one parameter family of 
solutions labeled by central field, and a given frequency 
is found for each value of the label as shown in uM ■ 
Excited solutions can be also constructed by allowing (f> 
to vanish at a given number of points, but always de- 
manding the satisfaction of the boundary conditions. Up 
to here the construction in spherical symmetry. Once 
we account with these data: i) we interpolated the wave 
function of the spherical data in the xz-grid and ii) re- 
solved the Poisson equation J5J , then we have initial data 
for ground state configurations in our axially symmetric 
domain. 

The system (|.':>l H is invariant under a clever scaling 
property given by 



{t,r,U,<f>} -> {A-^A-^A 2 ?/, A 2 0} (7) 
{p,M,K,W} -> {A 4 p, AM, A 3 X, A 3 IU} (8) 

where U is the gravitational potential, (f> is the spatial 
part of the wave function, p is the density of probability, 
M is the integral of p, K and W are the expectation value 
of the kinetic and gravitational energy respectively, and 



C. Testing the code with ground state 
configurations 

The steps followed in the construction of the solu- 
tions to equations llll'l) can be summarized as follows: 

1) choose a type of initial data for the t = slice if)(x, 0), 

2) populate the xz grid with those data, 3) solve equa- 
tion J2J) and get a gravitational potential, 4) using such 
potential leap the system using a At time slice fur- 
ther, 5) use the obtained wave function to solve again 
(0), 6) then repeat the loop 4, 5 either until the physics 
starts going wrong (physical quantities lose convergence, 
dissipative effects show up) or the cpu-time cannot be 
afforded. 

We test this code using systems whose properties we al- 
ready know from their construction, these are the ground 
state equilibrium configurations constructed in the pre- 
vious subsection; they show a particular property: in the 
continuum limit the wave function oscillates with con- 
stant frequency u>, which implies the density of probabil- 
ity p — \if)\ 2 and therefore the gravitational potential U 
remain time-independent. Unfortunately we do not ac- 
count with infinite resolution and therefore we are solv- 
ing the discretized versions of equations ill 121 ; instead we 
solve the second order finite differencing approximation 
of those equations, which is nothing but a truncated ex- 
pansion of the solution of the functions involved. Let us 
describe the situation with an example: assume we start 
with a ground state configuration with -0(0,0) = 1 and 
therefore the central density p(0,t) — p(0,0) = 1 at all 
times; first, as we are solving a truncated system of equa- 
tions the time-independence of the central density cannot 
be satisfied in a strict fashion, instead we can at most de- 
mand the central density to converge to 1 at all times; 
second, due to the truncation error of the finite differ- 
encing we are perturbing the system all the times and 
therefore the system should behave as an stable config- 
uration that is perturbed and thus should oscillate with 
the modes obtained from a perturbation theory analysis. 
Therefore we have a two-fold trap to verify whether or 
not the code is solving the physical system whose prop- 
erties we know beforehand. 

With all this in mind we evolved such a configuration 
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with tjj(0, 0) = 1 with our code and the results are as 
follows. The Fourier transform in Fig. Q]of the central 
density reveals that the main frequency of oscillation is 
7 = 0.046 which coincides with the result predicted by 
the first order radial perturbation theory developed in 
0; this indicates that despite the non spherical nature 
of our grid, the perturbation due to the discretization 
of the physical domain is spherical. On the right hand 
side of Fig. ^ we show the second order convergence 
of the central density to one, which indicates that our 
approximations work as they should when we refine the 
grid and approach the continuum limit. 




FIG. 1: Left: we show the Fourier transform of the cen- 
tral density of the configuration after the evolution has been 
performed; the main peak shows up at 7 = 0.046, which 
coincides with the result found using the perturbation the- 
ory in spherical symmetry from Right: the second or- 
der convergence of the central density to the value one is 
shown using two different resolutions. The runs were carried 
out on a x £ [0,20], z £ [—20,20] domain with resolutions 
Ax = Az = 0.1, 0.2. A Cauchy type convergence test decides 
whether or not we have convergence: the fact that the low 
resolution run shows a central density four times bigger than 
the one with the double resolution with respect to the value 
one, indicates the second order convergence to one. 
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FIG. 2: In this plot we show the meaning of virialized con- 
figurations. That is, due to discretization errors when solving 
the SP equations, there is an intrinsic error in the calcula- 
tions. We say the system is virialized only in the continuum 
limit. This is the reason why we show here the second order 
convergence also for the virial relation. Because this exam- 
ple corresponds to an ground state equilibrium configuration, 
which we know is virialized, we can be confident that not only 
the evolution code, but also the diagnostics tools work fine. 

In Fig. [3 we verify that ground state equilibrium con- 



figurations are virialized in the continuum limit, where 
the relation 2K + W — is satisfied with second order 
convergence. The quantities K and W are calculated as 
follows: 



K = --J i;*V 2 ^d 3 x (9) 

W = ^ ip*U^d 3 x (10) 

where the integrations are performed over the numerical 
domain. 



III. NON-SPHERICAL PERTURBATIONS 

It is still possible to use the discretization error to per- 
turb a ground state equilibrium configuration in a non- 
spherical way, for instance, using different resolutions in 
the x and z directions (see 0] for the use of such trick 
in relativistic boson stars). However this time we choose 
to fully-truly perturb the system with the addition of 
a shell of particles as done in 0, but this time using 
non-spherical shells like in [2(J. We start with a spher- 
ically symmetric ground state equilibrium configuration 
and add up a contribution proportional to a given spher- 
ical harmonic. Thus, we propose initial data given by 

V>(X, t) = tp ground + H (11) 

where dtp = e -*V* 2 +i/ 2 (^Lo a ' r z°)> where the co- 
efficients ai = Aie~ l -^ x2+z2 ~ r "^^' j2 are gaussian like 
shells added to the real part of the wave function of the 
ground state configuration; on the other hand, the factor 

e -it r yV+y contributes with a non-zero initial speed of 
the perturbation; notice that the up labels of the spheri- 
cal harmonics are all zero because otherwise they cannot 
be defined within an axially symmetric grid. Thus the 
complete set of parameters characterizing the perturba- 
tion are: ro, Ai, a and k r . 

We practiced several combinations of coefficients ai, 
including spherical perturbations containing only the Y§ 
contribution and found similar results. We only present 
the one with the parameters: ro = 10.0, Aq = Ai = 0, 
A2 = 0.02, a — 2.0, k r — —2.0 carried out on a grid 
x S [0, 20], z £ [-20, 20] with resolution Ax = Az = 0.2. 
The perturbation shows a quadrupolar contribution. In 
Fig. [3] we show the density of probability p at initial 
time. The mass of the shell is 0.8% that of the ground 
state equilibrium configuration. 

In Fig. Q]we show the evolution of such system. What 
can be seen is that the system relaxes and virializes 
around a spherical ground state equilibrium configura- 
tion. The original equilibrium configuration is recovered, 
which was verified in the following terms: the mass M 
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FIG. 3: The initial density profile of the perturbed equi- 
librium configuration. The central blob corresponds to the 
equilibrium configuration, whose magnitude is one at the ori- 
gin. The perturbation thus consists of two blobs coming with 
speed k r from the poles and a belt of particles over the equa- 
torial plane. 




FIG. 4: Left: evolution of the total energy; it can be observed 
that initially the whole system appears overwarmed due to 
the presence of the perturbation and its dynamical state; as 
the time runs, in a rather short time, the system recovers the 
total energy of the original configuration. Right: the quantity 
2K+W is monitored and it relaxes and converges to zero with 
second order. 



and the total energy E = K + W approach the values 
of the equilibrium configuration, the relation 2K + W 
converges to zero in the continuum limit after a short 
time, the system has a non-spherical initial shape and 
after a while the ellipticity converges to zero in the con- 
tinuum limit. This reasons indicate that these ground 
state equilibrium configurations are stable against non- 
spherical perturbations that involve the introduction of 
a quadrupolar shell of particles. 



A. Requirements for a Non-spherical Collapse 

The method followed to verify that non-spherical ini- 
tial configurations evolve toward a spherically symmetric 
ground state is as follows. 1) Given the evolution of a 
non-spherical initial profile ip(x,t) is carried out, we ob- 
tain via a Fourier Transform of a physical quantity, e.g. 
the central density pf , the frequency 7 of the fundamen- 
tal mode of oscillation of the system; We assume that 7 
corresponds to the characteristic frequency of oscillation 
of a linearly perturbed ground state ipperturbed toward 
which ip is evolving to. This is because the intrinsic per- 
turbation of the system associated to the discretization of 
the equations, we do not expect that, even for large evolu- 
tion times, tp matches exactly a stationary ground state. 
2) Once we estimate / (as done for the construction of 
Fig. Q we can calculate the rescaling parameter A that 
relates ipperturbed and the ground state with -0(x, 0) = 1. 
Using the definition of frequency and the rescaling rela- 
tion for t given by J7|) the rescaling parameter is calcu- 
lated as A = (J//) 1 / 2 . 3) The rescaling parameter A lets 
us calculate the non-hat quantities in (JSJ for the ground 
state "ipperturbed- Because we are not in the continuum 
limit we can at most demand that the physical quantities 
for ip(x, t), such as its density pf and mass M^, converge 
to those of a ground state configuration. 4) However, the 
information related to the convergence to an equilibrium 
configuration is not enough and we also verify that the 
fate of ip is a spherical configuration. We define the ellip- 
ticity of the system as the integrated difference between 
pf = p(0, z, t) and pf — p(x, 0, t) measured from the cen- 
ter of mass of the configuration. We observe for all our 
non-spherical initial configurations that after a transient 
period their ellipticity of the system relaxes and becomes 
zero in the continuum limit. 5) Finally, for all the evolved 
initial configurations we verify that the virialization con- 
dition is satisfied as shown in the perturbation case in 
Sec. EH 

B. Bigger perturbations 

In order to illustrate this we show the evolution of a 
ground state tp gr ound plus a considerable non-spherical 
density contribution which can be already considered to 
be a non-spherical initial profile. We evolved initial con- 
figurations of the form 



IV. NON-SPHERICAL COLLAPSE 

In this section we essay a step forward in the obvi- 
ous direction, that is, the collapse of non-spherical initial 
profiles for the SP system of equations. In fact we show 
that spherically symmetric ground states, are late time 
attractors for initial configurations which are not spher- 
ically symmetric. 



V>(x, t) = ipground + £ 2 =0 a/Y; (12) 

where a; = Ai exp[— r 2 /af] is a gaussian with r = 
V x 2 + z 2 centered at the origin, width 07 and real am- 
plitude Ai. Several runs were made for different A[ and 
07, and here we present only a representative one. In all 
cases the result is that the whole system evolves toward 
a ground state with a mass considerably bigger than that 
of tpground alone. The way in which the system evolves to 



a ground state is through the gravitational cooling pro- 
cess 0, which is powered by the ejection of scalar 
field/ 

This attractor behavior has been shown for spherically 
symmetric configurations in 13] and is shown here for 
the first time for the case of non-spherical initial profiles. 
Here we present the results for the initial configuration 
for which tpground (0, 0) = 2, I = 1, Ai = 1 and cr; = 2. 
The mass added to the system is the order ~ 3% that 
of the original ground state configuration, which implies 
the following result: some of the added particles joined 
the ground state configuration and evolved toward an- 
other rescalled ground state configuration different from 
the original one, simply because it is not a small pertur- 
bation this time. The evolution was carried out on a grid 
x e [0, 12], z e [-12, 12] with resolutions Ax = Az = 0.1 
and Ax — Az = 0.15. For this configuration we found 
the fundamental frequency at late time to be / = 0.0983, 
which gives a rescaling parameter of A = 1.462. The cen- 
tral density and mass of the ground state ip\ associated 
with this specific A are M\ = 3.014 and p\ — 4.563. In 
Fig. we show the evolution of the mass M$ versus the 
central density pf of the state tp for the two different res- 
olutions previously specified. The solid line is the branch 
of all the ground states constructed as indicated in sec- 
tionlrfland the cross symbol corresponds to the state with 
the parameters above. 

Convergence of and pf quantities to the star in 
the plot is shown also in Fig. [SJ The convergence is 
second order and the stars indicate the configuration we 
expect the system will relax onto for each resolution. A 
convergence analysis would reveal second order approach 
toward the configuration marked with the star. In Fig. 
EJ we show the ellipticity of the system; we can see that 
the system relaxes and becomes spherical. Also in this 
figui 
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FIG. 5: Evolution of the mass of the initial configuration 
versus the central density specified in the text. The resolu- 
tions used are Ax = Az — 0.1 (dotted) and Ax = Az — 0.15 
(dashed). The solid line represents the branch of spherical 
ground states, considering our calculations are second order 
convergent, from these two runs with the respective resolu- 
tions, we infer that the configuration in the cuntinuum limit 
is that marked with a cross. 
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FIG. 6: Evolution of the Ellipticity and the expression 2K + 
W. It is shown that after a while the initial axisymmetric 
configuration evolves toward a spherical one as the Ellipticity 
goes to zero. On the other hand as 2K + W oscillates around 
zero with decreasing amplitude we conclude that the system 
tends to a virialized state. 



C. Non-Spherical initial profiles 

What is in turn is to investigate the evolution of ax- 
isymmetric initial configurations. We choose these initial 
data to have the form 



1&(x, t) = Sf =0 A ; exp[-r 2 /a 2 ]lf (13) 

where r is as before, ai is the width of a gaussian and 
A\ a real amplitude. Several runs were made for differ- 
ent Ai and ai here we present the results for the initial 
configuration with Aq — 9.0, A\ = A-i = A3 — 1.0 and 
(Jo = 01 = 02 = 03 = 1.5. The evolution was carried 
out on a grid x £ [0,12], z G [—12,12] with resolution 
Ax = Az = 0.1 and Ax = Az = 0.15. For this con- 
figuration we found a fundamental frequency / = 0.128 
that implies a rescaling parameter A = 1.669. The cen- 
tral density and mass of the ground state <p\ associated 
with this specific A are M\ = 3.441 and p\ — 7.75 re- 
spectively, a point which is marked with a cross in Fig. 
[7| In such Figure we show the evolution of the mass 
versus the central density pf of the state -0 for the 
two grid resolutions previously specified. The solid line 
is the branch of all the ground states constructed as indi- 
cated in section m The stars are the configurations the 
system is approaching to for the different resolutions; be- 
cause we know such systems are only approximately an 
equilibrium configuration we practiced a Richardson ex- 
trapolation calculation and found that in terms of mass 
and central density, in the continuum limit -assuming the 
second order convergence we have in our algorithms- the 
cross is the configuration these results converge to. Thus 
in Figs. [7] and [H] we show that the system evolves toward 
a spherical and virialized configuration. 
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dxz=0.15 
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FIG. 7: Evolution of the central density and mass of the ini- 
tial configuration specified in the text. We used once again 
resolutions Ax — Az = 0.1 (dashed) and Ax = Az = 0.15 
(dotted). The solid line represents the branch of spherical 
ground states. As in the previous subsection, the stars corre- 
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FIG. 8: Temporal evolution of the Ellipticity and the ex- 
pression 2K + W of the initial configuration specified in the 
text. It is shown that after a while the initial axisymmetric 
configuration evolves toward a spherical one as the Ellipticity 
goes to zero. On the other hand as 2K + W oscillates around 
zero with decreasing amplitude we conclude that the system 
tends to a virialized state. 



with rather general axially symmetric shells of particles, 
which were proportional to the first axisymmetric spher- 
ical harmonics. Moreover, these shells were endowed 
with an initial speed toward the equilibrium configura- 
tion. The main result is that ground state configurations 
are stable against this type of perturbations. 

Finally, we have evolved axisymmetric initial profiles, 
and showed that spherically symmetric ground state 
equilibrium configurations play the role of late-time at- 
tractors. This is a generalization of the same result when 
the initial profiles are spherically symmetric ^3[- We 
have thus shown that the final state of a non-spherical 
self-gravitating scalar field fluctuation is a ground state 
spherically symmetric solution. Thus we showed that 
such configurations are late-time attractors not only for 
initial spherically symmetric configurations but also for 
quite general initial axisymmetric shapes. 

In the context of the scalar field dark matter model 
we have quite a new result: the collapse of overdensities 
tolerates an initial non-spherical contribution to the ini- 
tial profiles, and moreover, initially axisymmetric profiles 
tend toward a spherical ground state. 

About the process of virialization, with the results at 
hand we are unable to establish whether the gravitational 
cooling (powered by the emission of scalar field) is the 
only responsible for the relaxation of the system. The 
cases involving a non-zero quadrupolar contribution to 
the gravitational potential might allow one to formulate 
the study of another channel, that of the emission of grav- 
itational waves, which would be the subject of a further 
investigation. 

Definitely a major application of the tools shown in 
this paper is the demonstration of the stability properties 
of rotating solutions to the SP system of equations con- 
structed in the past 0, . It remains unclear whether 
excited state solutions representing spinning configura- 
tions arc stable, and whether such solutions could con- 
tribute to the SFDM possibilities in terms of quick virial- 
ization of collapsed structures and other dynamical prop- 
erties, for instance, the non-cupsy density profile might 
demand a considerably important spherical contribution 
in the case the wave function has odd parity. 



CONCLUSIONS 



We have presented a new code designed to solve the SP 
system with axial symmetry in cylindrical coordinates. 
We have shown it passes the necessary testbeds of sta- 
bility and consistency with expected results related to 
spherically symmetric ground state equilibrium configu- 
rations. 

We perturbed ground state equilibrium configurations 
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